scripts/Empirical example HADS/hads.R
In joakimwallmark/PolyOptimalIRT:
library(TestGardener)
library(mirt)
load(file="data/hads_dataList.rda")
load(file="data/hads_parList.rda")
# ----------- read in data -----------
# OUTDATED, SEE SPECIFIC DEPRESSION AND ANXIETY FILES
titlestr <- "hads-7-anxiety-screen"
U <- scan("data/hads.txt","o")
N <- length(U) # Number of examinees
Umat <- as.integer(unlist(stringr::str_split(U,"")))
n <- length(Umat)/N # Number of items
U <- matrix(Umat,N,n,byrow=TRUE)
# drop cases containing missing data
keepindex <- U[,14] < 4
N = sum(keepindex)
# N = 807
U <- U[keepindex,]
# gives names for mirt
colnames(U) <- paste("I", seq(1, ncol(U), 1), sep = "")
hadsGrm = mirt(data = U, model=1, itemtype = 'graded', SE=T)
hadsGrm2 = mirt(data = U, model=2, itemtype = 'graded', SE=T)
summary(hadsGrm2)
anova(hadsGrm, hadsGrm2) # p = 0 indicates better fit
#save(U, hadsGrm, file="data/hads_fittedmodel_mirt.RData")
for(i in 1:ncol(U)){
ItemPlot <- itemfit(hadsGrm,
group.bins=15,
empirical.plot = i,
empirical.CI = .95,
method = 'ML')
print(ItemPlot)
}
itemfit(hadsGrm)
# convert U from score format to index format by adding 1
U <- U + 1
# define key: NULL for scales
key <- NULL
# set scrtype, which is FALSE of the item is a rating scale
scrtype <- logical(n)
# define noption: 4 for 1:7 ,and 5 for 8 and 9 which containing missing values
noption <- matrix(4,n,1)
for (i in 1:n) noption[i] <- length(unique(U[,i]))
# dimension of space containing the curves
Wdim = sum(noption)
# summary count for each option
for (j in 1:n)
{
print(paste("Item: ",j,sep=""))
for (m in 1:noption[j])
{
print(paste("Option: ",m," N= ",sum(U[,j]==m),sep=""))
}
}
# --------- Define the option score values for each item ---------
ScoreList <- list()
for (item in 1:14){
ScoreList[[item]]<- c(0:3)
}
# compute the maximum sum score value
scrmax <- 0
for (j in 1:N) {
sumscrj = sum(U[j,]-1)
if (scrmax < sumscrj) scrmax <- sumscrj
}
scrmin <- 0
scrrng <- c(scrmin,scrmax)
# Set up the list vector containing the specifications of the data
itemVec <- c("Tense or wound-up", "Still enjoy things", "Frightened of something awful",
"Can laugh", "Worrying thoughts", "Cheerful", "Relaxed", "Slowed down",
"Frightened feeling", "Lost interest in appearance", "Restless",
"Look forward with enjoyment", "Panic", "Enjoy media")
optVec <-c('Not at all', 'Several days', 'More than half the days', 'Nearly every day')
optLab <- vector("list",n)
optLab[[ 1]] <- c("Most of the time", "A lot of the time", "From time to time", "Not at all")
optLab[[ 2]] <- c("Definitely as much", "Not quite so much", "Only a little", "Hardly at all")
optLab[[ 3]] <- c("Very definitely and quite badly", "Yes, but not too badly",
"A little, but it doesn't worry me", "Not at all")
optLab[[ 4]] <- c("As much as I always could", "Not quite so much now",
"Definitely not so much now", "Not at all")
optLab[[ 5]] <- c("A great deal of the time", "A lot of the time", "From time to time, but not too often ",
"Only occasionally")
optLab[[ 6]] <- c("Not at all", "Not often", "Sometimes", "Most of the time")
optLab[[ 7]] <- c("Definitely", "Usually", "Not often", "Not at all")
optLab[[ 8]] <- c("Nearly all the time", "Very often", "Sometimes", "Not at all")
optLab[[ 9]] <- c("Not at all", "Occasionally", "Quite often", "Very often")
optLab[[10]] <- c("Definitely", "I don't take as much care as I should ", "I may not take quite as much care",
"I take just as much care as ever ")
optLab[[11]] <- c("Very much indeed", "Quite a lot", "Not very much", "Not at all")
optLab[[12]] <- c("As much as I ever did", "Rather less than I used to",
"Definitely less than I used to", "Hardly at all")
optLab[[13]] <- c("Very often indeed", "Quite often", "Not very often", "Not at all")
optLab[[14]] <- c("Often", "Sometimes", "Not often", "Very seldom")
optList <- list(itemLab=itemVec, optLab=optLab, optScr=ScoreList)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr, NumBasis = 5, nbin=12)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr)
# save(hads_dataList, file="data/hads_dataList.rda")
# plot initial density
# ttllab <- paste(titlestr,": sum score", sep="")
# png(filename="hads_InitDensity.png")
# scoreDensity(hads_dataList$scrvec, hads_dataList$scrrng, ttlstr=ttllab)
# dev.off()
# --------- Set initial values that are required in the later analysis ---------
theta <- hads_dataList$percntrnk
thetaQnt <- hads_dataList$thetaQnt
# Bin the data, and smooth the binned data
chartList <- hads_dataList$chartList
WfdResult <- Wbinsmth(theta, hads_dataList, thetaQnt, chartList)
WfdList <- WfdResult$WfdList
binctr <- WfdResult$aves
Qvec <- c(5,25,50,75,95)
Wbinsmth.plot(binctr, Qvec, WfdList, hads_dataList, Wrng=c(0,3.5))
# --------------- Optimal scoring: cycle of smoothing/theta estimation ------------
# Set number of cycles and the cell array to containing the parameter
ncycle <- 20
# ----------------------------------------------------------------------------
# Proceed through the cycles
# ----------------------------------------------------------------------------
AnalyzeResult <- Analyze(theta, thetaQnt, hads_dataList, ncycle=ncycle, itdisp=FALSE)
AnalyzeResult.def <- AnalyzeResult
AnalyzeResult.bas5.bin12 <- AnalyzeResult
save(hads_dataList, AnalyzeResult, file = "data/hads_fittedmodel.RData")
parList <- AnalyzeResult$parList
MeanHvec <- AnalyzeResult$meanHvec
# ----------------------------------------------------------------------------
# Plot meanHsave and choose cycle for plotting
# ----------------------------------------------------------------------------
# cycleno <- 1:ncycle
# par(mfrow=c(1,1))
# hads_fittedmodel.RDatapng(filename="hads_Optimization.png")
# plot(cycleno,MeanHvec, type="b", lwd=2, xlab="Cycle Number")
# dev.off()
# select cycle for plotting
icycle <- 20
parListi <- parList[[icycle]]
WfdList <- parListi$WfdList
Qvec <- parListi$Qvec
binctr <- parListi$binctr
theta <- parListi$theta
arclength <- parListi$arclength
alfine <- parListi$alfine
# ----------------------------------------------------------------------------
# Plot surprisal curves for each test question
# ----------------------------------------------------------------------------
# plot both the probability and surprisal curves along with data points
hads_dataList$titlestr = "hads"
Wbinsmth.plot(binctr, Qvec, WfdList, hads_dataList, Wrng=c(0,4), saveplot=F)
# ----------------------------------------------------------------------------
# Plot density of theta
# ----------------------------------------------------------------------------
ttllab <- paste(titlestr,": percent rank", sep="")
edges <- c(0,100)
theta_in <- theta[0 < theta & theta < 100]
png(filename="hads_thetadensity.png")
scoreDensity(theta_in, edges, ttlstr=ttllab)
dev.off()
# ----------------------------------------------------------------------------
# Compute expected test scores for all examinees
# Plot expected test scores and expected test score over mesh
# ----------------------------------------------------------------------------
mu <- testscore(theta, WfdList, optList)
ttllab <- paste(titlestr,": expected score", sep="")
png(filename="hads_mudensity.png")
scoreDensity(mu, hads_dataList$scrrng, ttlstr=ttllab)
dev.off()
# compute expected score for each value in the fine mesh of theta values
indfine <- seq(0,100,len=101)
mufine <- testscore(indfine, WfdList, optList)
png(filename="hads_muplot.png")
mu.plot(mufine, hads_dataList$scrrng, titlestr=hads_dataList$titlestr)
dev.off()
# ----------------------------------------------------------------------------
# Compute arc length over a fine mesh of theta values and plot
# ----------------------------------------------------------------------------
# print length of the test info curve
print(paste("Arc length =", round(arclength,2)))
# plot arc length over fine mesh
png(filename="hads_arclenplot.png")
ArcLength.plot(arclength, alfine, titlestr=hads_dataList$titlestr)
dev.off()
# ----------------------------------------------------------------------------
# Display test effort curve projected into its first two principal components
# ----------------------------------------------------------------------------
# nharm=2
png(filename="hads_infoplot.png")
Result <- Wpca.plot(arclength, WfdList, Wdim, titlestr=hads_dataList$titlestr)
dev.off()
# nharm=3
Wpca.plot(arclength, WfdList, hads_dataList$Wdim, 3,
dodge = 1.005,titlestr=titlestr)
# ----------------------------------------------------------------------------
# Display sensitivity curves
# ----------------------------------------------------------------------------
# This code needs to put in a legend and indication of right answer if
# scoring is multiple choice
Sensitivity.plot(WfdList, Qvec, hads_dataList, titlestr=hads_dataList$titlestr,
width=c(-0.35,0.35), saveplot=T)
# ----------------------------------------------------------------------------
# Display power curves
# ----------------------------------------------------------------------------
Result = Power.plot(WfdList, Qvec, hads_dataList, saveplot=T)
# ----------------------------------------------------------------------------
# Display H, DH and D2H curves for selected examinee who has two minimal
# points, but only one of them has a strong positive 2nd derivative.
# ----------------------------------------------------------------------------
Hfuns.plot(theta, WfdList, hads_dataList$U, plotindex=8)
# ----------------------------------------------------------------------------
# Simulate data in order to show the root-mean-squared error,
# sampling error, and bias in estimates of theta and m(theta)
# ----------------------------------------------------------------------------
simList <- dataSimulation(hads_dataList, parListi, nsample=500)
dataSimulation.plot(simList)
# save.image(file="/Users/jimramsay/Documents/R/Umanitoba/hads/hads.rda")
joakimwallmark/PolyOptimalIRT documentation built on Dec. 21, 2021, 1:16 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(TestGardener)
library(mirt)
load(file="data/hads_dataList.rda")
load(file="data/hads_parList.rda")
# ----------- read in data -----------
# OUTDATED, SEE SPECIFIC DEPRESSION AND ANXIETY FILES
titlestr <- "hads-7-anxiety-screen"
U <- scan("data/hads.txt","o")
N <- length(U) # Number of examinees
Umat <- as.integer(unlist(stringr::str_split(U,"")))
n <- length(Umat)/N # Number of items
U <- matrix(Umat,N,n,byrow=TRUE)
# drop cases containing missing data
keepindex <- U[,14] < 4
N = sum(keepindex)
# N = 807
U <- U[keepindex,]
# gives names for mirt
colnames(U) <- paste("I", seq(1, ncol(U), 1), sep = "")
hadsGrm = mirt(data = U, model=1, itemtype = 'graded', SE=T)
hadsGrm2 = mirt(data = U, model=2, itemtype = 'graded', SE=T)
summary(hadsGrm2)
anova(hadsGrm, hadsGrm2) # p = 0 indicates better fit
#save(U, hadsGrm, file="data/hads_fittedmodel_mirt.RData")
for(i in 1:ncol(U)){
ItemPlot <- itemfit(hadsGrm,
group.bins=15,
empirical.plot = i,
empirical.CI = .95,
method = 'ML')
print(ItemPlot)
}
itemfit(hadsGrm)
# convert U from score format to index format by adding 1
U <- U + 1
# define key: NULL for scales
key <- NULL
# set scrtype, which is FALSE of the item is a rating scale
scrtype <- logical(n)
# define noption: 4 for 1:7 ,and 5 for 8 and 9 which containing missing values
noption <- matrix(4,n,1)
for (i in 1:n) noption[i] <- length(unique(U[,i]))
# dimension of space containing the curves
Wdim = sum(noption)
# summary count for each option
for (j in 1:n)
{
print(paste("Item: ",j,sep=""))
for (m in 1:noption[j])
{
print(paste("Option: ",m," N= ",sum(U[,j]==m),sep=""))
}
}
# --------- Define the option score values for each item ---------
ScoreList <- list()
for (item in 1:14){
ScoreList[[item]]<- c(0:3)
}
# compute the maximum sum score value
scrmax <- 0
for (j in 1:N) {
sumscrj = sum(U[j,]-1)
if (scrmax < sumscrj) scrmax <- sumscrj
}
scrmin <- 0
scrrng <- c(scrmin,scrmax)
# Set up the list vector containing the specifications of the data
itemVec <- c("Tense or wound-up", "Still enjoy things", "Frightened of something awful",
"Can laugh", "Worrying thoughts", "Cheerful", "Relaxed", "Slowed down",
"Frightened feeling", "Lost interest in appearance", "Restless",
"Look forward with enjoyment", "Panic", "Enjoy media")
optVec <-c('Not at all', 'Several days', 'More than half the days', 'Nearly every day')
optLab <- vector("list",n)
optLab[[ 1]] <- c("Most of the time", "A lot of the time", "From time to time", "Not at all")
optLab[[ 2]] <- c("Definitely as much", "Not quite so much", "Only a little", "Hardly at all")
optLab[[ 3]] <- c("Very definitely and quite badly", "Yes, but not too badly",
"A little, but it doesn't worry me", "Not at all")
optLab[[ 4]] <- c("As much as I always could", "Not quite so much now",
"Definitely not so much now", "Not at all")
optLab[[ 5]] <- c("A great deal of the time", "A lot of the time", "From time to time, but not too often ",
"Only occasionally")
optLab[[ 6]] <- c("Not at all", "Not often", "Sometimes", "Most of the time")
optLab[[ 7]] <- c("Definitely", "Usually", "Not often", "Not at all")
optLab[[ 8]] <- c("Nearly all the time", "Very often", "Sometimes", "Not at all")
optLab[[ 9]] <- c("Not at all", "Occasionally", "Quite often", "Very often")
optLab[[10]] <- c("Definitely", "I don't take as much care as I should ", "I may not take quite as much care",
"I take just as much care as ever ")
optLab[[11]] <- c("Very much indeed", "Quite a lot", "Not very much", "Not at all")
optLab[[12]] <- c("As much as I ever did", "Rather less than I used to",
"Definitely less than I used to", "Hardly at all")
optLab[[13]] <- c("Very often indeed", "Quite often", "Not very often", "Not at all")
optLab[[14]] <- c("Often", "Sometimes", "Not often", "Very seldom")
optList <- list(itemLab=itemVec, optLab=optLab, optScr=ScoreList)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr, NumBasis = 5, nbin=12)
hads_dataList <- make.dataList(U, key, optList, scrrng, titlestr)
# save(hads_dataList, file="data/hads_dataList.rda")
# plot initial density
# ttllab <- paste(titlestr,": sum score", sep="")
# png(filename="hads_InitDensity.png")
# scoreDensity(hads_dataList$scrvec, hads_dataList$scrrng, ttlstr=ttllab)
# dev.off()
# --------- Set initial values that are required in the later analysis ---------
theta <- hads_dataList$percntrnk
thetaQnt <- hads_dataList$thetaQnt
# Bin the data, and smooth the binned data
chartList <- hads_dataList$chartList
WfdResult <- Wbinsmth(theta, hads_dataList, thetaQnt, chartList)
WfdList <- WfdResult$WfdList
binctr <- WfdResult$aves
Qvec <- c(5,25,50,75,95)
Wbinsmth.plot(binctr, Qvec, WfdList, hads_dataList, Wrng=c(0,3.5))
# --------------- Optimal scoring: cycle of smoothing/theta estimation ------------
# Set number of cycles and the cell array to containing the parameter
ncycle <- 20
# ----------------------------------------------------------------------------
# Proceed through the cycles
# ----------------------------------------------------------------------------
AnalyzeResult <- Analyze(theta, thetaQnt, hads_dataList, ncycle=ncycle, itdisp=FALSE)
AnalyzeResult.def <- AnalyzeResult
AnalyzeResult.bas5.bin12 <- AnalyzeResult
save(hads_dataList, AnalyzeResult, file = "data/hads_fittedmodel.RData")
parList <- AnalyzeResult$parList
MeanHvec <- AnalyzeResult$meanHvec
# ----------------------------------------------------------------------------
# Plot meanHsave and choose cycle for plotting
# ----------------------------------------------------------------------------
# cycleno <- 1:ncycle
# par(mfrow=c(1,1))
# hads_fittedmodel.RDatapng(filename="hads_Optimization.png")
# plot(cycleno,MeanHvec, type="b", lwd=2, xlab="Cycle Number")
# dev.off()
# select cycle for plotting
icycle <- 20
parListi <- parList[[icycle]]
WfdList <- parListi$WfdList
Qvec <- parListi$Qvec
binctr <- parListi$binctr
theta <- parListi$theta
arclength <- parListi$arclength
alfine <- parListi$alfine
# ----------------------------------------------------------------------------
# Plot surprisal curves for each test question
# ----------------------------------------------------------------------------
# plot both the probability and surprisal curves along with data points
hads_dataList$titlestr = "hads"
Wbinsmth.plot(binctr, Qvec, WfdList, hads_dataList, Wrng=c(0,4), saveplot=F)
# ----------------------------------------------------------------------------
# Plot density of theta
# ----------------------------------------------------------------------------
ttllab <- paste(titlestr,": percent rank", sep="")
edges <- c(0,100)
theta_in <- theta[0 < theta & theta < 100]
png(filename="hads_thetadensity.png")
scoreDensity(theta_in, edges, ttlstr=ttllab)
dev.off()
# ----------------------------------------------------------------------------
# Compute expected test scores for all examinees
# Plot expected test scores and expected test score over mesh
# ----------------------------------------------------------------------------
mu <- testscore(theta, WfdList, optList)
ttllab <- paste(titlestr,": expected score", sep="")
png(filename="hads_mudensity.png")
scoreDensity(mu, hads_dataList$scrrng, ttlstr=ttllab)
dev.off()
# compute expected score for each value in the fine mesh of theta values
indfine <- seq(0,100,len=101)
mufine <- testscore(indfine, WfdList, optList)
png(filename="hads_muplot.png")
mu.plot(mufine, hads_dataList$scrrng, titlestr=hads_dataList$titlestr)
dev.off()
# ----------------------------------------------------------------------------
# Compute arc length over a fine mesh of theta values and plot
# ----------------------------------------------------------------------------
# print length of the test info curve
print(paste("Arc length =", round(arclength,2)))
# plot arc length over fine mesh
png(filename="hads_arclenplot.png")
ArcLength.plot(arclength, alfine, titlestr=hads_dataList$titlestr)
dev.off()
# ----------------------------------------------------------------------------
# Display test effort curve projected into its first two principal components
# ----------------------------------------------------------------------------
# nharm=2
png(filename="hads_infoplot.png")
Result <- Wpca.plot(arclength, WfdList, Wdim, titlestr=hads_dataList$titlestr)
dev.off()
# nharm=3
Wpca.plot(arclength, WfdList, hads_dataList$Wdim, 3,
dodge = 1.005,titlestr=titlestr)
# ----------------------------------------------------------------------------
# Display sensitivity curves
# ----------------------------------------------------------------------------
# This code needs to put in a legend and indication of right answer if
# scoring is multiple choice
Sensitivity.plot(WfdList, Qvec, hads_dataList, titlestr=hads_dataList$titlestr,
width=c(-0.35,0.35), saveplot=T)
# ----------------------------------------------------------------------------
# Display power curves
# ----------------------------------------------------------------------------
Result = Power.plot(WfdList, Qvec, hads_dataList, saveplot=T)
# ----------------------------------------------------------------------------
# Display H, DH and D2H curves for selected examinee who has two minimal
# points, but only one of them has a strong positive 2nd derivative.
# ----------------------------------------------------------------------------
Hfuns.plot(theta, WfdList, hads_dataList$U, plotindex=8)
# ----------------------------------------------------------------------------
# Simulate data in order to show the root-mean-squared error,
# sampling error, and bias in estimates of theta and m(theta)
# ----------------------------------------------------------------------------
simList <- dataSimulation(hads_dataList, parListi, nsample=500)
dataSimulation.plot(simList)
# save.image(file="/Users/jimramsay/Documents/R/Umanitoba/hads/hads.rda")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.